Major questions:
Firstly, we’re going to install and load our packages. We will use the pacman Package Management Tool which allows us to install and load subsequent packages in a condensed and efficient way.
# Install and load required packages
pacman::p_load(dplyr, ggmap, ggspatial, ggplot2, ggpubr, janitor, lme4, lmerTest, lwgeom, MuMIn, readxl, sf, sp, tidyr, tidyverse, viridis)Secondly, we’ll read in our data and explore it.
Let’s view the data.
Now, we’ll wrangle some of our data to get it ready for analyses and plotting.
data <- read_excel(path="data.xlsx", sheet="data") %>%
clean_names() %>%
mutate(date = as.Date(date)) %>%
mutate(bait_take = as.numeric(bait_take)) %>%
mutate(phase = as.factor(phase)) %>%
drop_na(bait_take) %>%
mutate(bait_yes = ifelse(bait_take == 1, "Yes", "No"))Ellie, explain why you’re going to exclude Fridays here
The weather data were acquired from the Australian Bureau of Meteorology weather station 70351.
# Subset to actual instances of bait take by a vertebrate
data <- data %>%
subset(species != "Lizard") %>%
mutate(fox_visited = ifelse(fox_visit_1 != 0, "Yes",
ifelse(fox_visit_2 !=0, "Yes",
ifelse(fox_visit_3 !=0, "Yes",
ifelse(fox_visit_4 !=0, "Yes", "No"))))) %>%
mutate(fox_visits = ifelse(fox_visit_4 != 0, 4,
ifelse(fox_visit_3 !=0, 3,
ifelse(fox_visit_2 !=0, 2,
ifelse(fox_visit_1 !=0, 1, 0))))) %>%
mutate(raven_visited = ifelse(raven_visit_1 != 0, "Yes",
ifelse(raven_visit_2 !=0, "Yes",
ifelse(raven_visit_3 !=0, "Yes", "No")))) %>%
mutate(raven_visits = ifelse(raven_visit_3 !=0, 3,
ifelse(raven_visit_2 !=0, 2,
ifelse(raven_visit_1 !=0, 1, 0))))
write.csv(data, "data processed.csv")Bait take by treatment
# Read in processed data
data <- read.csv("data processed.csv") %>%
clean_names() %>%
mutate(date = as.Date(date))
# Plot by treatment type
bait_treat <- ggplot(data, aes(x = as.Date(date),
y = bait_take,
col=treatment,
fill=treatment)) +
geom_smooth() +
theme_minimal() +
theme(plot.title = element_text(hjust=0.5)) +
ggtitle("All species") +
xlab("Date") +
ylab("Bait take") +
labs(fill="Treatment", col="Treatment")
# Display the plot
print(bait_treat)Bait take by treatment and phase
# Read in processed data
data <- read.csv("data processed.csv") %>%
clean_names() %>%
mutate(date = as.Date(date))
# Plot by treatment type
bait_treat_phase <- ggplot(data, aes(x = as.Date(date),
y = bait_take,
col=treatment,
fill=treatment)) +
geom_smooth() +
theme_minimal() +
theme(plot.title = element_text(hjust=0.5)) +
ggtitle("All species") +
facet_wrap(~phase, scales="free") +
xlab("Date") +
ylab("Bait take") +
labs(fill="Treatment", col="Treatment")
# Display the plot
print(bait_treat_phase)Percentage bait taken by each species by phase
# Subset to only potential foxes
taken <- data %>%
subset(bait_take != 0) %>%
subset(species != "Not taken") %>%
group_by(species, treatment, phase) %>%
summarise(perc = n() / nrow(data) * 100)
# Plot by treatment type
taken_species_phase_hist <- ggplot(data=taken,
aes(x = species, y = perc,
col=treatment, fill=treatment)) +
geom_col() +
theme_minimal() +
theme(plot.title = element_text(hjust=0.5)) +
facet_wrap(~phase) +
ggtitle("All species") +
xlab("Species") +
ylab("Percentage of taken baits (%)") +
labs(fill="Treatment", col="Treatment")
# Display the plot
print(taken_species_phase_hist)Plot bait take by rainfall
# Plot by rainfall type
bait_rain <- ggplot(data, aes(x = rainfall_mm,
y = bait_take,
col=treatment,
fill=treatment)) +
geom_smooth() +
theme_minimal() +
theme(plot.title = element_text(hjust=0.5)) +
ggtitle("All species") +
xlab("Rainfall") +
ylab("Bait take") +
labs(fill="Treatment", col="Treatment")
# Display the plot
print(bait_rain)Percentage fox visits by treatment and phase
# Subset to only potential foxes
fox_stats <- data %>%
subset(species == "Fox") %>%
group_by(treatment, phase) %>%
summarise(perc = n() / nrow(data) * 100)
# Plot by treatment type
visits_fox_histogram <- ggplot(data=fox_stats,
aes(y = perc, x = treatment,
col=treatment, fill=treatment)) +
geom_col() +
theme_minimal() +
theme(plot.title = element_text(hjust=0.5)) +
facet_wrap(~phase) +
ggtitle("Potential fox visits (including unknowns)") +
xlab("Species") +
ylab("Percentage of fox visits (%)") +
labs(fill="Treatment", col="Treatment")
# Display the plot
print(visits_fox_histogram)# Subset to only potential foxes
foxes_only <- data %>%
subset(species != "Raven") %>%
subset(species != "Unknown")
# Plot by treatment type
bait_treat_phase <- ggplot(data=foxes_only, aes(x = as.Date(date),
y = bait_take,
col=treatment,
fill=treatment)) +
geom_smooth() +
theme_minimal() +
theme(plot.title = element_text(hjust=0.5)) +
ggtitle("Foxes vs no take (excluding unknowns)") +
facet_wrap(~phase, scales="free") +
xlab("Date") +
ylab("Bait take") +
labs(fill="Treatment", col="Treatment")
# Display the plot
print(bait_treat_phase)Percentage fox visits by treatment and phase
# Subset to only potential foxes
raven_stats <- data %>%
subset(species == "Raven") %>%
group_by(treatment, phase) %>%
summarise(perc = n() / nrow(data) * 100)
# Plot by treatment type
visits_raven_hist <- ggplot(data=raven_stats,
aes(y = perc, x = treatment,
col=treatment, fill=treatment)) +
geom_col() +
theme_minimal() +
theme(plot.title = element_text(hjust=0.5)) +
facet_wrap(~phase) +
ggtitle("Potential raven visits (including unknowns)") +
xlab("Species") +
ylab("Percentage of raven visits (%)") +
labs(fill="Treatment", col="Treatment")
# Display the plot
print(visits_raven_hist)# Subset to only potential foxes
ravens_only <- data %>%
subset(species != "Fox") %>%
subset(species != "Unknown")
# Plot by treatment type
bait_raven_treat_phase <- ggplot(data=ravens_only, aes(x = as.Date(date),
y = bait_take,
col=treatment,
fill=treatment)) +
geom_smooth() +
theme_minimal() +
theme(plot.title = element_text(hjust=0.5)) +
ggtitle("Ravens vs no take (excluding unknowns)") +
facet_wrap(~phase, scales="free") +
xlab("Date") +
ylab("Bait take") +
labs(fill="Treatment", col="Treatment")
# Display the plot
print(bait_raven_treat_phase)Bait take by fox visits
# Plot bait take by species visits
bait_visits <- ggplot(data, mapping=aes(x = as.Date(date))) +
geom_smooth(data, mapping=aes(y = as.numeric(bait_take)),
color = "red", fill = "red") +
geom_smooth(data, mapping=aes(y = as.numeric(fox_visits)),
color = "darkorange", fill = "darkorange") +
geom_smooth(data, mapping=aes(y = as.numeric(raven_visits)),
color = "black", fill = "black") +
theme_minimal() +
facet_wrap(~treatment) +
scale_y_continuous(name = "Bait take (%)", limits = c(0, 1),
sec.axis = sec_axis(~., name = "Fox (orange) and raven (black) visits",
breaks = c(0, 1))) +
xlab("Date")
# Display the plot
print(bait_visits)Build generalised linear models for our research questions.
##
## Call:
## glm(formula = bait_take ~ treatment, data = data)
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.66038 0.02223 29.705 < 2e-16 ***
## treatmentTreatment -0.10124 0.03164 -3.199 0.00142 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for gaussian family taken to be 0.2357498)
##
## Null deviance: 224.02 on 941 degrees of freedom
## Residual deviance: 221.60 on 940 degrees of freedom
## AIC: 1316.1
##
## Number of Fisher Scoring iterations: 2
# Model for bait take by treatment and rainfall
summary(glm(bait_take ~ treatment + rainfall_mm, data=data))##
## Call:
## glm(formula = bait_take ~ treatment + rainfall_mm, data = data)
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.638871 0.023237 27.493 < 2e-16 ***
## treatmentTreatment -0.100640 0.031505 -3.194 0.00145 **
## rainfall_mm 0.006263 0.002060 3.040 0.00243 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for gaussian family taken to be 0.2337005)
##
## Null deviance: 224.02 on 941 degrees of freedom
## Residual deviance: 219.44 on 939 degrees of freedom
## AIC: 1308.9
##
## Number of Fisher Scoring iterations: 2
##
## Call:
## glm(formula = bait_take ~ treatment + phase, data = data)
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.69484 0.05679 12.234 < 2e-16 ***
## treatmentTreatment -0.10146 0.03165 -3.205 0.00139 **
## phase -0.01787 0.02710 -0.659 0.50980
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for gaussian family taken to be 0.2358916)
##
## Null deviance: 224.02 on 941 degrees of freedom
## Residual deviance: 221.50 on 939 degrees of freedom
## AIC: 1317.7
##
## Number of Fisher Scoring iterations: 2
##
## Call:
## glm(formula = bait_take ~ treatment + site, data = data)
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.901039 0.148702 6.059 2.01e-09 ***
## treatmentTreatment -0.102424 0.065102 -1.573 0.116011
## siteF01 0.046330 0.173125 0.268 0.789062
## siteF02 0.201385 0.170912 1.178 0.238991
## siteF03 -0.014532 0.174168 -0.083 0.933523
## siteF04 -0.157455 0.168687 -0.933 0.350858
## siteF05 -0.134257 0.170471 -0.788 0.431161
## siteF06 -0.051039 0.171986 -0.297 0.766718
## siteF07 0.098961 0.170001 0.582 0.560632
## siteF08 -0.003736 0.173070 -0.022 0.982781
## siteF09 0.161306 0.167225 0.965 0.335005
## siteF10 -0.025888 0.171773 -0.151 0.880239
## siteF11 -0.003736 0.173070 -0.022 0.982781
## siteF12 -0.032747 0.169541 -0.193 0.846886
## siteF13 -0.466256 0.169131 -2.757 0.005956 **
## siteF14 -0.560520 0.172711 -3.245 0.001216 **
## siteF15 -0.585677 0.170318 -3.439 0.000611 ***
## siteF16 -0.374846 0.168687 -2.222 0.026524 *
## siteF17 -0.310130 0.170001 -1.824 0.068444 .
## siteF18 -0.480433 0.171773 -2.797 0.005270 **
## siteF19 -0.621452 0.171156 -3.631 0.000298 ***
## siteF20 -0.571342 0.171773 -3.326 0.000916 ***
## siteF21 -0.651039 0.171986 -3.785 0.000164 ***
## siteF22 -0.610447 0.170471 -3.581 0.000361 ***
## siteF23 -0.466256 0.169131 -2.757 0.005956 **
## siteF24 -0.578201 0.169541 -3.410 0.000678 ***
## siteF25 -0.453736 0.173070 -2.622 0.008898 **
## siteF26 -0.434979 0.171773 -2.532 0.011502 *
## siteF27 -0.466256 0.169131 -2.757 0.005956 **
## siteF28 0.201385 0.170912 1.178 0.238991
## siteF29 0.196729 0.171156 1.149 0.250693
## siteF30 0.106147 0.172711 0.615 0.538981
## siteF31 -0.232064 0.172070 -1.349 0.177788
## siteF32 0.058162 0.169541 0.343 0.731635
## siteF34 -0.084329 0.172711 -0.488 0.625480
## siteF35 0.008052 0.170001 0.047 0.962232
## siteF40 -0.666907 0.171156 -3.896 0.000105 ***
## siteF48 -0.751374 0.174168 -4.314 1.78e-05 ***
## siteF52 0.104719 0.167760 0.624 0.532646
## siteF53 0.103414 0.168687 0.613 0.539996
## siteF54 -0.537745 0.170912 -3.146 0.001708 **
## siteF55 -0.091515 0.170949 -0.535 0.592552
## siteF56 0.008601 0.170471 0.050 0.959774
## siteF57 -0.526493 0.167171 -3.149 0.001690 **
## siteF58 -0.810130 0.170001 -4.765 2.20e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for gaussian family taken to be 0.1493372)
##
## Null deviance: 223.41 on 937 degrees of freedom
## Residual deviance: 133.36 on 893 degrees of freedom
## (4 observations deleted due to missingness)
## AIC: 924.16
##
## Number of Fisher Scoring iterations: 2
# Do fox visits change with treatment or phase?
summary(glm(fox_visits ~ treatment + phase +
rainfall_mm + maximum_temperature_c, data=data))##
## Call:
## glm(formula = fox_visits ~ treatment + phase + rainfall_mm +
## maximum_temperature_c, data = data)
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.181054 0.226700 -0.799 0.4247
## treatmentTreatment 0.022190 0.050005 0.444 0.6573
## phase 0.036098 0.053793 0.671 0.5024
## rainfall_mm 0.007035 0.003418 2.058 0.0398 *
## maximum_temperature_c 0.014803 0.009938 1.490 0.1367
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for gaussian family taken to be 0.5886676)
##
## Null deviance: 556.46 on 941 degrees of freedom
## Residual deviance: 551.58 on 937 degrees of freedom
## AIC: 2181.1
##
## Number of Fisher Scoring iterations: 2
# Model for visits against treatment and phase
summary(glm(raven_visits ~ treatment + phase +
rainfall_mm + maximum_temperature_c, data=data))##
## Call:
## glm(formula = raven_visits ~ treatment + phase + rainfall_mm +
## maximum_temperature_c, data = data)
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.120731 0.101268 1.192 0.233
## treatmentTreatment 0.016744 0.022338 0.750 0.454
## phase -0.002641 0.024030 -0.110 0.913
## rainfall_mm -0.002171 0.001527 -1.422 0.155
## maximum_temperature_c -0.001225 0.004439 -0.276 0.783
##
## (Dispersion parameter for gaussian family taken to be 0.1174665)
##
## Null deviance: 110.37 on 941 degrees of freedom
## Residual deviance: 110.07 on 937 degrees of freedom
## AIC: 662.88
##
## Number of Fisher Scoring iterations: 2
Here we map our responses of bait_take,
fox_visits, and raven_visits across
Ginninderry Conservation Corridor.
# Use Google API to fetch a base map
# ggmap::register_google(key="[add API key here]", write=TRUE)
map <- get_map(location=c(lon=148.9811, lat=-35.2350),
zoom=15, source="google", maptype="satellite", crop=FALSE)# Read in shapefile
gis <- st_read("shapefiles/bait_stations_egg_cta.shp") %>%
fortify() %>%
clean_names() %>%
rename(site=name)## Reading layer `bait_stations_egg_cta' from data source
## `X:\Honours\R\shapefiles\bait_stations_egg_cta.shp' using driver `ESRI Shapefile'
## Simple feature collection with 43 features and 5 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: 148.9673 ymin: -35.24309 xmax: 148.9927 ymax: -35.22524
## Geodetic CRS: GDA2020
# Map the bait stations
gin_map <- ggmap(map) +
geom_point(data=gis, size=2, aes(x=long, y=lat, col=treatment)) +
theme_minimal() +
xlab("") + ylab("") + labs(col="Treatment")
# Display the plot
print(gin_map)# Combine the bait take df with the GIS df
data <- left_join(data, gis, by="site")
# Generalised linear model
mod <- glm(bait_take ~ lat + long + elevation, data=data)
summary(mod)##
## Call:
## glm(formula = bait_take ~ lat + long + elevation, data = data)
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -4.841e+03 7.170e+02 -6.753 2.56e-11 ***
## lat -1.071e+01 3.202e+00 -3.345 0.000857 ***
## long 2.997e+01 4.978e+00 6.020 2.51e-09 ***
## elevation -7.231e-05 6.421e-04 -0.113 0.910369
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for gaussian family taken to be 0.1887321)
##
## Null deviance: 222.12 on 930 degrees of freedom
## Residual deviance: 174.95 on 927 degrees of freedom
## (11 observations deleted due to missingness)
## AIC: 1095.7
##
## Number of Fisher Scoring iterations: 2
# Map the bait stations
mod_map <- ggmap(map) +
geom_point(data=data, shape=21, size=2, stroke=1, col="white", alpha=0.2,
aes(x=long, y=lat, fill=bait_take)) +
theme(axis.text.x=element_blank(),
axis.text.y=element_blank(),
axis.ticks.x=element_blank(),
axis.ticks.y=element_blank(),
plot.title=element_text(hjust=0.5),
strip.background=element_rect(fill="#363842"),
strip.text=element_text(colour="white")) +
facet_wrap(~phase) +
scale_fill_viridis_c() +
xlab("") + ylab("") + ggtitle("Phase") + labs(fill="Bait take")
# Display the plot
print(mod_map)# Map the fox visitation
mod_map <- ggmap(map) +
geom_point(data=data, shape=21, size=2, stroke=1, col="white", alpha=0.2,
aes(x=long, y=lat, fill=fox_visits)) +
theme(axis.text.x=element_blank(),
axis.text.y=element_blank(),
axis.ticks.x=element_blank(),
axis.ticks.y=element_blank(),
plot.title=element_text(hjust=0.5),
strip.background=element_rect(fill="#363842"),
strip.text=element_text(colour="white")) +
facet_wrap(~phase) +
scale_fill_viridis_c() +
xlab("") + ylab("") + ggtitle("Phase") + labs(fill="Fox visits")
# Display the plot
print(mod_map)# Map the fox visitation
mod_map <- ggmap(map) +
geom_point(data=data, shape=21, size=2, stroke=1, col="white", alpha=0.2,
aes(x=long, y=lat, fill=raven_visits)) +
theme(axis.text.x=element_blank(),
axis.text.y=element_blank(),
axis.ticks.x=element_blank(),
axis.ticks.y=element_blank(),
plot.title=element_text(hjust=0.5),
strip.background=element_rect(fill="#363842"),
strip.text=element_text(colour="white")) +
facet_wrap(~phase) +
scale_fill_viridis_c() +
xlab("") + ylab("") + ggtitle("Phase") + labs(fill="Raven visits")
# Display the plot
print(mod_map)Here we read in, project, and map the bait stations across Ginninderry Conservation Corridor (GCC) using the following shapefiles:
What we need to do for each layer is check its inherent Coordinate
Reference System (CRS, using st_transform() to find the
last 4-digit number of its output), remind the layer of that number
(using st_set_crs()), and then reproject it to the CRS we
wish to work with (using st_transform("EPSG:4326")).
# Read in bait station points
stations <- st_read("shapefiles/bait_stations_egg_cta.shp") %>%
as.data.frame() %>%
fortify() %>%
clean_names()## Reading layer `bait_stations_egg_cta' from data source
## `X:\Honours\R\shapefiles\bait_stations_egg_cta.shp' using driver `ESRI Shapefile'
## Simple feature collection with 43 features and 5 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: 148.9673 ymin: -35.24309 xmax: 148.9927 ymax: -35.22524
## Geodetic CRS: GDA2020
# Tranform to Spatial Points Dataframe
stations_spdf <- SpatialPointsDataFrame(
data.frame(stations$long, stations$lat),
stations, proj4string=CRS("EPSG:4326")) %>%
st_as_sf()## Reading layer `gcc_boundary_polygon' from data source
## `X:\Honours\R\shapefiles\gcc_boundary_polygon.shp' using driver `ESRI Shapefile'
## Simple feature collection with 1 feature and 1 field
## Geometry type: POLYGON
## Dimension: XY
## Bounding box: xmin: 148.9616 ymin: -35.24409 xmax: 148.9934 ymax: -35.21299
## Geodetic CRS: WGS 84
# Remind the layer of its CRS then reproject to our desired CRS
gcc_bou <- gcc_bou %>%
# Set projection using sf
st_set_crs(st_crs(gcc_bou)) %>%
st_transform("EPSG:4326")
# Calculate distance of bait stations to nearest edge of polygon
dist_bou <- st_geometry(obj = gcc_bou) %>%
st_cast(to = 'MULTILINESTRING') %>%
st_distance(y=stations_spdf) %>%
as.data.frame() %>%
# Transpose the df from columns to rows
t()
# Join distances to bait stations df
stations_spdf <- cbind(stations_spdf, dist_bou)## Reading layer `i5516_roads' from data source
## `X:\Honours\R\shapefiles\i5516_roads.shp' using driver `ESRI Shapefile'
## Simple feature collection with 2775 features and 17 fields
## Geometry type: MULTILINESTRING
## Dimension: XY
## Bounding box: xmin: 148.5 ymin: -36 xmax: 150 ymax: -35
## Geodetic CRS: GDA94
# Remind the layer of its CRS then reproject to our desired CRS
gcc_roads <- gcc_roads %>%
# Set projection using sf
st_set_crs(st_crs(gcc_roads)) %>%
st_transform("EPSG:4326")
# Convert linestrings to a single multilinestring
gcc_roads <- st_union(gcc_roads)
# Calculate distance of bait stations to nearest edge of polygon
dist_roads <- st_geometry(obj = gcc_roads) %>%
st_cast(to = 'MULTILINESTRING') %>%
st_distance(y=stations_spdf) %>%
as.data.frame() %>%
# Transpose the df from columns to rows
t()
# Join distances to bait stations df
stations_spdf <- cbind(stations_spdf, dist_roads)## Reading layer `gcc_construction_area' from data source
## `X:\Honours\R\shapefiles\gcc_construction_area.shp' using driver `ESRI Shapefile'
## Simple feature collection with 1 feature and 1 field
## Geometry type: POLYGON
## Dimension: XY
## Bounding box: xmin: 148.9725 ymin: -35.22814 xmax: 148.9903 ymax: -35.21997
## Geodetic CRS: GDA2020
# Remind the layer of its CRS then reproject to our desired CRS
construct <- construct %>%
# Change EPSG here to match the last number from st_crs() output
st_set_crs(st_crs(construct)) %>%
st_transform("EPSG:4326")
# Calculate distance of bait stations to nearest edge of polygon
dist_con <- st_geometry(obj = construct) %>%
st_cast(to = 'MULTILINESTRING') %>%
st_distance(y=stations_spdf) %>%
as.data.frame() %>%
# Transpose the df from columns to rows
t()
# Join distances to bait stations df
stations_spdf <- cbind(stations_spdf, dist_con)## Reading layer `i5516_nativevegetationareas' from data source
## `X:\Honours\R\shapefiles\i5516_nativevegetationareas.shp' using driver `ESRI Shapefile'
## Simple feature collection with 313 features and 13 fields
## Geometry type: POLYGON
## Dimension: XY
## Bounding box: xmin: 148.5 ymin: -36 xmax: 150 ymax: -35
## Geodetic CRS: GDA94
# Remind the layer of its CRS then reproject to our desired CRS
gcc_veg <- gcc_veg %>%
# Change EPSG here to match the last number from st_crs() output
st_set_crs(st_crs(gcc_veg)) %>%
st_transform("EPSG:4326")
# Check for and fix invalid geometries
gcc_veg <- st_make_valid(gcc_veg)
# Convert multiple features to a single feature
gcc_veg <- st_union(gcc_veg)
# Calculate distance of bait stations to nearest edge of polygon
dist_veg <- st_geometry(obj = gcc_veg) %>%
st_cast(to = 'MULTILINESTRING') %>%
st_distance(y=stations_spdf) %>%
as.data.frame() %>%
# Transpose the df from columns to rows
t()
# Join distances to bait stations df
stations_spdf <- cbind(stations_spdf, dist_veg)## Reading layer `i5516_pondageareas' from data source
## `X:\Honours\R\shapefiles\i5516_pondageareas.shp' using driver `ESRI Shapefile'
## Simple feature collection with 6 features and 13 fields
## Geometry type: POLYGON
## Dimension: XY
## Bounding box: xmin: 149.0126 ymin: -35.31815 xmax: 149.5838 ymax: -35.05082
## Geodetic CRS: GDA94
# Remind the layer of its CRS then reproject to our desired CRS
gcc_pond <- gcc_pond %>%
# Change EPSG here to match the last number from st_crs() output
st_set_crs(st_crs(gcc_pond)) %>%
st_transform("EPSG:4326")
# Convert multiple features to a single feature
gcc_pond <- st_union(gcc_pond)
# Calculate distance of bait stations to nearest edge of polygon
dist_pond <- st_geometry(obj = gcc_pond) %>%
st_cast(to = 'MULTILINESTRING') %>%
st_distance(y=stations_spdf) %>%
as.data.frame() %>%
# Transpose the df from columns to rows
t()
# Join distances to bait stations df
stations_spdf <- cbind(stations_spdf, dist_pond)## Reading layer `i5516_watercourselines' from data source
## `X:\Honours\R\shapefiles\i5516_watercourselines.shp' using driver `ESRI Shapefile'
## Simple feature collection with 3627 features and 15 fields
## Geometry type: MULTILINESTRING
## Dimension: XY
## Bounding box: xmin: 148.5 ymin: -36.00002 xmax: 150 ymax: -35
## Geodetic CRS: GDA94
# Remind the layer of its CRS then reproject to our desired CRS
gcc_wcl <- gcc_wcl %>%
# Change EPSG here to match the last number from st_crs() output
st_set_crs("EPSG:4283") %>%
st_transform("EPSG:4326")
# Convert linestrings to a single multilinestring
gcc_wcl <- st_union(gcc_wcl)
# Calculate distance of bait stations to nearest edge of polygon
dist_wcl <- st_geometry(obj = gcc_wcl) %>%
st_cast(to = 'MULTILINESTRING') %>%
st_distance(y=stations_spdf) %>%
as.data.frame() %>%
# Transpose the df from columns to rows
t()
# Join distances to bait stations df
stations_spdf <- cbind(stations_spdf, dist_wcl)Here we model the relationship between our response variables (i.e.,
bait_take, fox_vists, and
raven_visits) and the distance between the bait station and
nearby spatial features (i.e., boundary [dist_bou], roads
[dist_roads], construction site [dist_con],
vegetation strip [dist_veg], pondage
[dist_pond], watercourses [dist_wcl]).
# Subset to only variables in the subsequent models
data_mod <- data %>%
select(bait_take, treatment, fox_visits, raven_visits,
dist_bou, dist_roads, dist_con, dist_veg, dist_pond, dist_wcl) %>%
na.omit()
# Fit the model with glm
mod <- glm(bait_take ~ treatment + dist_bou + dist_roads + dist_con +
dist_veg + dist_pond + dist_wcl,
data = data_mod, family=binomial, na.action = "na.fail")
# Display summary statistics
summary(mod)##
## Call:
## glm(formula = bait_take ~ treatment + dist_bou + dist_roads +
## dist_con + dist_veg + dist_pond + dist_wcl, family = binomial,
## data = data_mod, na.action = "na.fail")
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.7788067 7.0046783 -0.397 0.691583
## treatmentTreatment -0.9047092 0.1729551 -5.231 1.69e-07 ***
## dist_bou -0.0033458 0.0008050 -4.156 3.23e-05 ***
## dist_roads 0.0030916 0.0005053 6.119 9.42e-10 ***
## dist_con 0.0032323 0.0009439 3.424 0.000616 ***
## dist_veg 0.0016249 0.0007593 2.140 0.032349 *
## dist_pond -0.0002584 0.0012654 -0.204 0.838209
## dist_wcl 0.0003089 0.0006154 0.502 0.615655
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 1247.77 on 930 degrees of freedom
## Residual deviance: 966.61 on 923 degrees of freedom
## AIC: 982.61
##
## Number of Fisher Scoring iterations: 4
# Since dist_pond and dist_wcl was non-significant, remove from model
mod <- glm(bait_take ~ treatment + dist_bou + dist_roads + dist_con + dist_veg,
data = data_mod, family=binomial, na.action = "na.fail")
# Display summary statistics
summary(mod)##
## Call:
## glm(formula = bait_take ~ treatment + dist_bou + dist_roads +
## dist_con + dist_veg, family = binomial, data = data_mod,
## na.action = "na.fail")
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -4.2406271 0.5364397 -7.905 2.68e-15 ***
## treatmentTreatment -0.8817769 0.1671349 -5.276 1.32e-07 ***
## dist_bou -0.0031113 0.0006613 -4.704 2.55e-06 ***
## dist_roads 0.0030243 0.0004669 6.478 9.29e-11 ***
## dist_con 0.0031826 0.0003490 9.118 < 2e-16 ***
## dist_veg 0.0018967 0.0002137 8.877 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 1247.8 on 930 degrees of freedom
## Residual deviance: 966.9 on 925 degrees of freedom
## AIC: 978.9
##
## Number of Fisher Scoring iterations: 4
## Global model call: glm(formula = bait_take ~ treatment + dist_bou + dist_roads +
## dist_con + dist_veg, family = binomial, data = data_mod,
## na.action = "na.fail")
## ---
## Model selection table
## (Int) dst_bou dst_con dst_rds dst_veg trt df logLik AICc delta
## [1,] -4.241 -0.003111 0.003183 0.003024 0.001897 + 6 -483.448 979 0
## weight
## [1,] 1
## Models ranked by AICc(x)
Plot bait take against dist_bou, dist_con,
dist_roads, dist_veg.
# Plot the distances against bait take
dist_bou_plot <- ggplot(data = data_mod) +
geom_smooth(aes(x=dist_bou, y=bait_take),
col="purple", fill="purple") +
facet_wrap(~treatment) +
xlab("Distance to nearest boundary (m)") +
ylab("Bait take (proportion)") +
theme_minimal()
dist_con_plot <- ggplot(data = data_mod) +
geom_smooth(aes(x=dist_con, y=bait_take),
col="orange", fill="orange") +
facet_wrap(~treatment) +
xlab("Distance to nearest part of construction site (m)") +
ylab("Bait take (proportion)") +
theme_minimal()
dist_road_plot <- ggplot(data = data_mod) +
geom_smooth(aes(x=dist_roads, y=bait_take),
col="darkgrey", fill="darkgrey") +
facet_wrap(~treatment) +
xlab("Distance to nearest road (m)") +
ylab("Bait take (proportion)") +
theme_minimal()
dist_veg_plot <- ggplot(data = data_mod) +
geom_smooth(aes(x=dist_veg, y=bait_take),
col="green", fill="green") +
facet_wrap(~treatment) +
xlab("Distance to nearest vegetation patch (m)") +
ylab("Bait take (proportion)") +
theme_minimal()
dist_plots <- ggarrange(dist_bou_plot, dist_con_plot,
dist_road_plot, dist_veg_plot,
nrow=2, ncol=2)
# Display the plot
print(dist_plots)Next, repeat these models and plots for fox visits and raven visits.
Map these layers
# Map the site
gcc_map <- ggplot() +
geom_sf(data=gcc_roads, col="darkgrey") +
geom_sf(data=construct, fill="orange", alpha=0.3) +
geom_sf(data=gcc_veg, fill="green", alpha=0.3) +
geom_sf(data=gcc_pond, col="lightblue") +
geom_sf(data=gcc_wcl, col="blue") +
geom_sf(data=gcc_bou, fill="purple", alpha=0.3) +
geom_sf(data=stations_spdf, aes(col=treatment)) +
geom_text(data=stations_spdf, aes(x=long, y=lat, label=name), size=3, color="black") +
# Limit map to GCC extent
coord_sf(xlim = c(st_bbox(stations_spdf)[[1]]-0.005,
st_bbox(stations_spdf)[[3]]+0.005),
ylim = c(st_bbox(stations_spdf)[[2]]-0.005,
st_bbox(stations_spdf)[[4]]+0.005)) +
theme_minimal() +
scale_fill_manual(values=viridis(18, begin=0.9, end=0.1)) +
xlab("") + ylab("") + labs(col="Treatment")
# Display the plot
print(gcc_map)Session information